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We study models of correlated percolation where there are constraints on the occupation of 
sites that mimic force-balance, i.e. for a site to be stable requires occupied neighboring sites in 
all four compass directions in two dimensions. We prove rigorously that p c < 1 for the two- 
dimensional models studied. Numerical data indicate that the force-balance percolation transition 
is discontinuous with a growing crossover length, with perhaps the same form as the jamming 
percolation models, suggesting the same underlying mechanism driving the transition in both cases. 
In other words, force-balance percolation and jamming percolation may indeed belong to the same 
universality class. We find a lower bound for the correlation length in the connected phase and 
that the correlation function does not appear to be a power law at the transition. Finally, we study 
the dynamics of the culling procedure invoked to obtain the force-balance configurations and find a 
dynamical exponent similar to that found in sandpile models. 



I. INTRODUCTION 

Uncorrelated percolation and its associated geometric 
phase transition is arguably the most studied paradigm 
for a phase transition in a disordered system. Not only 
have physicists been able to nail down the universality 
class of percolation, more recently, mathematicians have 
been able to rigorously verify the universality class for 
at least one particular percolation model [1, 2]. Even 
more recently, this work has been extended to a sec- 
ond two-dimensional percolation model [3]. Other two- 
dimensional models are expected to follow. 

While uncorrelated percolation exhibits a continuous 
phase transition, there exists a new class of correlated 
percolation models that provably exhibit a discontinuous 
phase transition — a notable departure from the transi- 
tion in uncorrelated models [4-8]. What do we mean 
by correlated percolation? We mean that there are con- 
straints imposed on the occupation of sites such that cor- 
relations in the occupation arise regardless of whether or 
not there is a transition. 

One of the simplest class of models of correlated perco- 
lation is fc-core/bootstrap percolation [9-13]. It is defined 
as follows. Consider a regular lattice of coordination 
number Z max , and some integer k with 2 < k < Z max . 
Initially, sites are independently occupied with probabil- 
ity p. Then, all occupied sites with fewer than k neigh- 
boring occupied sites are eliminated. This decimation 
process is repeatedly applied to the surviving occupied 
sites, until all surviving sites (if any) have at least k 
surviving neighbors. The surviving sites are called the 
fc-core, and phases of the model are determined by the 
presence or absence of an infinite cluster of these sur- 
vivors. The /c-core percolation model has a number of 
physical realizations [14], including nonmagnetic impuri- 
ties in a magnetic system [9] and the glass transition via 
a kinctically-constraincd spin-flip model known as the 
Fredrickson- Andersen model [15, 16]. 

Another well-known model of correlated percolation is 
the Kob- Andersen model [17], a particle-conserving coun- 
terpart to the Fredrickson-Andersen model. In the Kob- 



Andcrscn model, a particle can hop if and only if there 
are at least m empty neighbors before and after a par- 
ticle hop. As the density of particles increases, it be- 
comes more difficult to hop, and the density of frozen 
particles increases, resulting in slower dynamics. Even- 
tually, the frozen particles percolate throughout the sys- 
tem, resulting in a glass transition. Based on numerical 
analysis of the percolation of frozen particles, it was ini- 
tially thought that the Kob- Andersen model exhibited a 
glass transition at a value of p c ss 0.84 in two dimen- 
sions. However, recent work by Toninclli, Biroli, and 
Fisher rigorously demonstrates that the thermodynamic 
p c is actually unity for the range of m relevant to the 
glass transition [18, 19]. 

Jamming percolation, a new class of two-dimensional 
correlated percolation models inspired by kinetically con- 
strained models, was introduced by Toninclli, Biroli, and 
Fisher [4, 5, 8]. These models consist of the following oc- 
cupation constraints: a site can remain occupied only if 
there exists at least one occupied site in set A and one oc- 
cupied site in set B, or one occupied site in set C and one 
occupied site in set D. All sets are disjoint from each other 
and contain two sites. See Fig. 1 for the sets in a jamming 
percolation model called the spiral model. Jamming per- 
colation models exhibit discontinuous phase transitions 
with a crossover length scale that diverges faster than 
any power law. It has been recently demonstrated that 
this behaviour prevails when more than four disjoint sets 
are introduced [20]. While there arc rigorous results for 
some models in this class, it is unknown how generic this 
transition is and whether or not there are other unusual, 
or atypical, phase transitions of correlated percolation. 

Here, we present a class of correlated percolation mod- 
els denoted as force-balance percolation models. Force- 
balance percolation was originally introduced in Ref. [11] 
as a toy model for the jamming transition in finite di- 
mensions. The jamming transition is a transition from a 
liquid-like to an amorphous solid-like state as some par- 
ticular control parameter, such as the packing density, is 
varied. Examples of potentially related jamming systems 
arc glass-forming liquids, colloidal suspensions, foams, 
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FIG. 1: Sets for the spiral model. 



emulsions, and granular matter [21]. Despite decades 
of study of the glass-forming liquids in particular, how- 
ever, it is even unclear whether these transitions are true 
thermodynamic transitions or merely examples of kinetic 
arrest [22]. 

There has been some recent activity focusing on a zero- 
temperature jamming transition in a system of repulsive 
soft spheres as the packing density is increased. Numeri- 
cal simulations by O'Hern, et al. considered repulsive soft 
particles in two and three dimensions [23, 24]. For small 
packing fraction cf>, the particles easily arrange themselves 
so as not to overlap with any other particle, and the to- 
tal potential energy thus vanishes. As <fi is increased, 
there is a particular value of (j) c (Point J) above which 
the particles can no longer avoid each other and the to- 
tal potential energy becomes nonzero. The system jams 
in that it develops nonzero static bulk and shear moduli 
above cf> c . The average coordination number (the aver- 
age number of overlapping neighbors per particle) jumps 
from Z = to Z = Z c at Point J, and then rises with 
increasing packing fraction <fi as Z — Z c ~ (</> — <p c )P , 
where j3 — 0.5. This behaviour was recently observed in 
a two-dimensional systems of glass beads [25]. Further- 
more, as <p approaches <j) c from above, the singular part 
of the shear modulus vanishes with the exponent 7 = 0.5; 
more recent simulations also find that there is a length 
scale that diverges with an exponent v = 0.25 [26]. So 
the transition at Point J appears to have characteristics 
of both first-order and second-order phase transitions: at 
the transition, Z is discontinuous, but there are nontriv- 
ial power laws. 

To understand this somewhat unusual transition, an 
analogy to fc-core/bootstrap percolation was made in 
Ref. [11]. The scalar aspect of the principle of local me- 
chanical stability, where particles need d + 1 contacts, 
maps to a requirement of fc occupied neighbors. Surpris- 
ing agreement was found between the mean-field fc-core 
percolation exponents and the repulsive soft sphere sim- 
ulations. We note that the two- and three- dimensional 
simulations observed the same exponents, suggesting a 
possible critical dimension of two since logarithmic cor- 
rections would be difficult to determine [27]. 



While there is agreement between the mean field fc-core 
exponents and the low-dimensional repulsive soft sphere 
exponents, fc-core in finite-dimensional spaces does not 
appear to have such agreement. Such systems seem to fall 
into one of two classes: the transition is either continuous 
and in the same universality class as normal percolation, 
or it does not occur until p c = 1. In the first class, sys- 
tems that allow finite clusters all exhibit continuous tran- 
sitions [28, 29]. In the second class, large voids are very 
likely to grow, and in the infinite system limit, with prob- 
ability one there will be at least one void that will grow 
to empty the entire system [30, 31]. This prevents fc-core 
percolation for any p < 1. Unfortunately, neither cat- 
egory of fc-core percolation describes the nontrivial dis- 
continuous transition observed in the finite-dimensional 
simulations and experiments of jamming. 

In an effort to try to capture the behaviour of jam- 
ming in finite dimensions, force-balance percolation was 
introduced in Ref. [11]. In this model, the fc-core con- 
straint is retained, but the vectorial constraints of the 
principle of local mechanical stability are mimicked by 
creating culling rules that take into account where the 
neighboring particles are located. Loosely speaking, if 
there is a neighboring contact to one side of the particle, 
there must be at least one neighboring contact on the 
other side of the particle to allow for force balance. 

We must emphasize the force-balance percolation is a 
model with no explicit forces. We look only at connec- 
tivity, in contrast to models such as rigidity percolation 
where repulsive and attractive forces are defined on, say, 
a lattice of springs [32-35]. However, the nature and pos- 
sible universality of the rigidity percolation transition is 
still up for debate, despite decades of study. Since force- 
balance percolation is a much simpler model there is ulti- 
mately a better chance of analyzing it beyond numerics. 

Here, we explore several versions of force-balance 
percolation, both rigorously and numerically, to begin 
to answer the following questions: 

(0) Is there a force-balance percolation transition? 

(1) If there is indeed a transition, what are its proper- 
ties? Is it continuous? 

(2) How generic is the transition among the various 
force-balance models? 

(3) Is there a link between force-balance percolation and 
jamming percolation? Are they in the same universality 
class? 

The paper is organized as follows. In section II we re- 
view the force-balance percolation model introduced in 
Ref. [11], and introduce two related models. We present 
in section III a rigorous proof that the thermodynamic p c 
is less than unity for at least two of these models. Ear- 
lier work on fc-core percolation misinterpreted numerical 
results, finding transitions in novel universality classes, 
with p c < 1 [30, 36, 37], when in fact the models studied 
had p c = 1 [31, 38, 39]; our proof renders our subse- 
quent interpretation of the numerical data presented in 



section IV on sounder footing. Finally, we close in sec- 
tion V with a summary of our findings and discuss their 
implications. 



II. MODELS 

For the first force-balance percolation model, we begin 
with a two-dimensional square lattice. Each site neigh- 
bors all sites except itself within a 5x5 square — each site 
therefore has 24 nearest neighbors. (For a 3x3 square, 
with the following rules, p c — 1). Since we are in two 
dimensions, we impose a 3-core constraint. However, we 
also impose the force-balance constraint, which is the fol- 
lowing: there must be at least one occupied neighbor in 
set A, which in turn calls for at least one occupied neigh- 
bor in set B, and there must be at least one occupied 
neighbor in set C, which in turn calls for at least one 
occupied neighbor in set D. The four sets A, B, C, and 
D, are defined in Fig. 2. The force-balance constraint 
can be succinctly stated as: (A and B) and (C and D), 
where each letter X is short for "at least one occupied site 
in set X". Note that the force-balance constraint is de- 
fined in such a way such that vertical and/or horizontal 
lines of occupied particles are, by themselves, not sta- 
ble. Fig. 3a demonstrates an allowed configuration and 
Fig. 3b demonstrates a prohibited configuration. 

To enforce the force-balance and fc-core constraints, we 
initially occupy sites on the lattice with independent oc- 
cupation probabilities p, and then repeatedly remove oc- 
cupied sites that violate either the fc-core or force-balance 
constraints, until all remaining occupied sites are stable. 
Note that p is the occupation density before culling, and 
gcnerically differs from the final occupation density. 
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FIG. 2: Force-balance model on the 2d square lattice with 24 
nearest neighbors. 



The model is abelian. In other words, the final con- 
figuration after the culling process is independent of the 
order in which sites are culled. It can be done in parallel, 
or in series, or some combination thereof. One can define 
another force-balance-like model that allows for horizon- 
tal and vertical lines of occupied neighbors. However, 
such a model would be nonabelian. Nonabliean models 
are less desirable, because in such models the final re- 
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FIG. 3: Black sites denote occupied sites, and red sites de- 
note unoccupied sites, (a) Allowed configuration, (b) For- 
bidden configuration. While the number of occupied nearest 
neighbors is greater than three, the force-balance condition is 
violated. 



suits depend on the order in which sites are culled. In 
our work, the culling procedure is merely an algorithm 
to achieve the force-balance percolation configuration. 

To determine whether or not the behaviour observed 
in the original force-balance model is generic, we define 
two additional abelian models. Our second model is de- 
fined on the square lattice, with 16 nearest neighbors, 
and quadrants as defined in Fig. 4. Again, k — 3 and the 
force-balance constraint is the same as above: (A and B) 
and (C and D). This model is obviously similar to the 
first one, though the ratio of k to Z max is different. 

The third model is a three-dimensional model with 26 
nearest neighbors and six regions, color-coded in Fig. 5. 
Each of the six regions consists of a 3 x 3 square on a face 
of the 3x3x3 cube centered on the site whose stability 
is being analyzed. For this model, we impose a 4-core 
constraint (since d = 3) , and the force-balance condition 
requires an ocupied site in each of the six regions. 

For comparison, along with numerical simulations of 
these three models, we also performed simulations on the 
spiral model. This model was introduced by Toninelli, 
et al. [7, 8], and they have proved that it undergoes a 
jamming transition, and obtained rigorous results about 
the properties of the transition [8] . 
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FIG. 4: A second force-balance model on the 2d square lattice 
with 16 nearest neighbors. 
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FIG. 5: Force-balance model on the 3d cubic lattice with 26 
nearest neighbors. A red site participates in only one of the 
six sets, a purple site in two of the six sets, and a blue site in 
three of the six sets. 



III. RIGOROUS RESULTS 



Proof that p c < 1 in the 2d force-balance 
models 



probability of the origin participating in an infinite chain 
of 3-clusters to the southwest is also nonzero. Looking 
only at the infinite chain of 3-clusters, it is straightfor- 
ward to check that every site in it is stable under the 
culling rules. To see this, observe that each 3-cluster in 
the infinite path has only four possible configurations of 
adjacent 3-clusters, shown in Fig. 8. For all four con- 
figurations, all sites in the central 3-cluster are stable 
under the k = 6, NE-SW stability condition. Therefore, 
Pc < iPc P ) 1 ^ < 1- This bound holds for k < 6 for the 
force-balance model, since the force-balance constraints 
are a subset of the constraints in the NE-SW model. 

A similar proof holds for the 16 nearest-neighbor force- 
balance model. For this model, one can construct 2- 
clusters to obtain a proof that for k < 4 we have 

Pc<( P ? P ) 1/2 - 

We note that a simple lower bound for p c for both mod- 
els is the p c of the corresponding uncorrelated percolation 
models (i.e. the p c of the model without culling). 



To prove that p c < 1 for the initial force-balance model 
on the square lattice, we will demonstrate that p c < 1 
for a more heavily constrained model with the same 24 
nearest neighbors. In particular, we require (1) k = 6, 
(2) at least one occupied neighbor in the 4-site region to 
the northeast, and (3) at least one occupied neighbor in 
the 4-site region to the southwest. See Fig. 6 for the two 
regions. Since any sites stable under these conditions 
are automatically stable in the force-balance model, if 
p c < 1 for the k — 6 NE-SW model, then also p c < 1 for 
the force-balance model. 



NE 



SW 




FIG. 7: Adjacent clusters of 3 sites in the NE-SW model. 



FIG. 6: k = 6 NE-SW model. 



Next, we prove that the origin has a non-zero probabil- 
ity of participating in an infinite cluster for the k — 6 NE- 
SW model. We divide the lattice into clusters of three 
sites, as shown in Fig. 7. Certain clusters are defined 
as adjacent, and connected by directed lines in Fig. 7. 
Sites not grouped in 3-clusters play no role in the proof, 
and can be ignored. Now, suppose that p > (p^' p ) 1 ^ 3 , 
where p® p is the critical probability for two-dimensional 
directed percolation. Then each 3-cluster is occupied 
with probability p 3 > p^ p '. The directed lines in Fig. 
6 are isomorphic to two-dimensional percolation. Thus, 
the probability of the origin participating in an infinite 
chain of 3-clusters to the northeast is nonzero, provided 
the 3-cluster at the origin is occupied. Similarly, the 
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FIG. 8: Local configurations of 3-clusters in the NE-SW 
model. 



Regarding the three-dimensional model, as was shown 
in Ref. [20], for three-dimensional versions of the jam- 
ming percolation models, percolating structures along 
one direction can cross without having sites in common. 
This prevents an obvious mapping to directed percola- 
tion. We have therefore not obtained a rigorous bound 
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for p c in the three-dimensional model. 



G' 



B. Nature of the transition? 



While it can be shown that p c < 1 for the two- 
dimensional models, there is no proof as to whether the 
transition is continuous or discontinuous. For the two- 
dimensional jamming percolation models, there exists a 
rigorous argument for a discontinuous transition assum- 
ing a conjecture regarding directed percolation [4, 8]. 
The argument relics on having two independent directed 
percolation processes arising from the disjoint pairs of 
sets, each containing two sites. The jamming percolation 
transition can then be shown to occur at the directed 
percolation transition point. In the force-balance models 
considered here, the sets are composed of more than two 
sites, and are not disjoint. The critical occupation prob- 
ability thus cannot be easily identified. This means that 
it is not obvious how to construct a rigorous argument 
along the lines of the jamming percolation models. 

However, we first argue that the alteration of some 
properties of the jamming percolation models to make 
them look more like force-balance percolation should not 
change the nature of the jamming percolation transition. 
For example, the property of having two sites per set 
in the jamming percolation models can be extended to 
having more sites per set by investigating other lattice 
models with more nearest neighbors, but belonging to the 
same universality class as directed percolation. Surely, 
there exist other directed percolation processes with more 
than two neighbors per site. 

On the other hand, the force-balance rules can be mod- 
ified to be more like the jamming percolation models, 
where the "and" between overlapping sets is equivalent 
to "or" between various pairs or triplets of smaller dis- 
joint sets. See Fig. 9. While the occupation of pairs of 
disjoint sets are the same as in the jamming percolation 
case (excluding the fc-core constraint where more than 
one site in the set must be occupied for the pairs), the 
occupation of triplets are more stringent. It remains to 
be seen whether or not a rigorous argument can be con- 
structed for force-balance percolation. 

Given these arguments, the transition for the force bal- 
ance percolation models may in fact behave similarly to 
the jamming percolation models and so we expect the 
transition in the force-balance case to be discontinuous 
as well. We will numerically test this hypothesis as well 
as others in the upcoming section. 



FIG. 9: The force-balance rules can be implemented with (A' 
and C) or (D' and B') or (E' and B' and C) or (F' and C 
and D') or (G' and A' and B') or (ET and A' and D') or (A' 
and G' and FT) or (B' and E' and G') or (C and E' and F') 
or (D' and F' and H') or (E' and F' and G' and FT). 



IV. NUMERICAL RESULTS 
A. 24 NN model 

1. Culling dynamics 

We begin our numerical simulations of the 24 NN 
(nearest neighbor) force-balance model by looking at the 
dynamics of the culling process. The lattice is size L by 
L. We surround it with a boundary of fully occupied 
sites, two layers thick, which are never culled. We de- 
note this as wired boundary conditions. Then, at each 
timestep, we simultaneously cull every unstable site. We 
repeat this until all remaining sites are stable and record 
the culling time. 

For any given system size, the mean culling time as a 
function of the initial occupation density peaks for some 
density — a sample curve for L = 2896 is shown in Fig- 
ure 10. The density at which the curve peaks provides 
one possible definition of the critical density. The size 
of the peak is denoted by M, the maximum number of 
average culls. M grows with L. We plot M vs. L on a 
log- log scale in Figure 11. The results are well fit by a 
power law, although there is visible overall curvature. A 
best fit for L > 32 gives M cx L a , with a = 1.226 ±0.027, 
which is suggestively close to 5/4. 

We note that 5/4 is the dynamical exponent for 
avalanches in the sandpile model [40]. Priezzhcv has 
given an analytical argument that the dynamics expo- 
nent for the sandpile model should be exactly 5/4 [41]. 
More generally, Pietronero, et al. [40] have given a renor- 
malization group argument that the avalanche exponent 
should be the same for a wide class of "sandpile-like" 
models [42]. They obtain an approximate exponent of 
1.253 w 5/4. "Sandpile-like" models are those in which 
energy is dissipated, so that instabilities propagate from 
site to site. A similar process may be responsible for the 
exponent of 5/4 seen in our models here: culling at a 
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site may be analogous to toppling at a critical point of a 
sandpile-like model, in the sense that the culling of one 
site triggers adjacent cullings that can propagate long 
distances. 
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FIG. 10: Mean culling time as a function of initial occupation 
density for the 24 NN model, L = 2896. 
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FIG. 11: Peak of the mean culling time as a function of L for 
the 24 NN model. The best fit has a slope of 1.226 ± 0.027. 
The same data is also plotted for the spiral model. 



and no finite clusters are allowed, the culling avalanches 
should become spatially long-ranged near p c = 1. This 
result is to be contrasted with fc-core cases with finite fc- 
core clusters. There, no power law distribution of culling 
avalanche sizes was found [43]. 

In the force-balance model, there are no finite force- 
balance clusters, so the culling avalanches should be spa- 
tially long-ranged near the transition. Whether or not 
there should be a broad distribution of sizes is not clear 
a priori. Given the sandpile-like behavior detected in 
the mean culling time required to obtain a force-balance 
cluster, one may expect to find a broad distribution as is 
found in sandpile models. However, if the force-balance 
percolation transition is discontinuous one would expect 
a well-defined avalanche size for those systems whose re- 
dundant sites have already been removed. 

Numerically, we calculate the probability of having a 
culling avalanche size s, P(s) in the presence of peri- 
odic boundary conditions. See Fig. 12. The probability 
is broad near and above the transition for intermediate 
avalanche sizes. On a log-log plot, the slope of P(s) for 
the broadly distributed intermediate-sized avalanches de- 
pends somewhat on p and on L. More careful study 
is needed to determine whether or not these avalanche 
sizes are consistent with the measured —1.253 for sand- 
pile models [42]. It is of note that for p = 0.45 and 
L = 512, for example, the slope is approximately —1.3. 

As opposed to quantitative analysis of the 
intermediate-sized avalanches, we are more inter- 
ested in first determining the qualitative nature of the 
curves. As opposed to a purely broad distribution, there 
is a prominent peak at the tail of the distribution. These 
well-defined avalanche sizes correspond the last occuring 
avalanches when the lattice ultimately empties out. 
These correspond to the marginal infinite cluster. The 
peak persists as p is decreased towards the p c indicated 
from the position of the maximum in the mean culling 
times, suggesting a discontinuous transition for the onset 
of the infinite cluster. This behaviour is retained in the 
larger systems with the peak becoming more separated 
from the broad part of the distribution. See Fig. 13. 



2. Force-balance avalanches 

Farrow, et al. [43] studied the dynamics of culling 
in fc-core percolation on various lattices by looking at 
avalanches. More precisely, the culling process is iterated 
until a stable fc-core configuration is obtained. Then, a 
random site in the stable fc-core is removed. It triggers 
the removal of other sites until the fc-core stabilizes once 
more. The number of sites removed during this removal 
is the culling avalanche size. The process is repeated 
until the lattice is completely empty In cases where 
p c = 1 and the stability rules allow no finite clusters, 
a power-law distribution of culling avalanche sizes was 
found. For these cases, the system goes from being unoc- 
cupied for p c < 1 , to being fully occupied at p c = 1. Since 
all sites must be eventually be removed below p c = 1 



3. Spanning cluster 

We next look at P s , the probability of spanning. We 
define this, for wired boundary conditions, as the proba- 
bility that the largest cluster connects either the top and 
bottom sides, or the right and left sides; for this defini- 
tion, sites are defined as connected if one is within the 
24 site neighborhood of the other (see Fig. 2). Fig. 14 
depicts P s for the 24 NN model with wired boundary 
conditions. 

Fig. 15 shows the same curves, but for periodic bound- 
ary conditions. For periodic boundary conditions, it is 
impossible to have any occupied sites without having a 
spanning cluster, so to check if we have a spanning clus- 
ter, we just need to check that at least one site is occu- 
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FIG. 14: The probability of spanning for wired boundary 
conditions in the 24 NN model. 



FIG. 12: Log-log plot of P(s) for L — 128 using periodic 
boundary conditions. 
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FIG. 15: The probability of spanning for periodic boundary 
conditions in the 24 NN model. Error bars are too small to 
be seen. 



FIG. 13: Log-log plot of P(s) for L = 128 and L = 512 for 
comparison. The L = 512 curves have been shifted downward 
by a factor of 0.1. 



pied. 

For periodic boundary conditions, the P a curves can 
be generated quickly by a simple modification of the 
Newman-Ziff algorithm [44]. For uncorrelated percola- 
tion, the Newman-Ziff algorithm passes once through ev- 
ery possible density, in increasing order, by adding oc- 
cupied sites in random order, and updating the connec- 
tivity with the Hoshen-Kopelman algorithm [45]. This 
allows for efficient generation of P s for every occupa- 
tion density. This method is not available to us for 
the force-balance model with wired boundary conditions, 
because there sites are added by increasing the density, 
but also removed by culling. This prevents a single pass 
through all densities from being done in time Oil? In L). 
(The Hoshen-Kopelmann algorithm does not allow clus- 
ter identifications to be rapidly updated after removal of 
sites.) However, with periodic boundary conditions in a 
force-balance model, we do not need to run the Hoshen- 



Kopelman algorithm. We can instead start from a fully 
occupied lattice, and remove sites until the culling condi- 
tion causes an empty lattice. This allows the calculation 
of P s at every density in time 0(L 2 ). 

From Figs. 14 and 15, we can define the critical prob- 
ability for a specified system size, p c (L), as the initial 
occupation density that gives P s — 1/2. This definition 
differs from that based on the peak of the mean culling 
time, given in subsection IV A 1. However, as shown in 
figure 16, the two definitions have the same qualitative 
dependence on L, and approach each other for large sys- 
tem sizes. 

Given p c (L) we can extract the thermodynamic critical 
point, p c (L — oo), if given the functional dependence 
of p c (L) on L. However, it is unclear what functional 
dependence to fit to. It was assumed in Ref. [11] that 
\p c (L = oo) — p c (L)\ oc i" 1 /", where v is an exponent 
characterizing the divergence of the length scale. With 
this functional form, v can be detected by varying p c (oo) 
until a log-log plot gives a straight line. However, such a 
procedure is not a check of the functional form, since we 
are allowed to choose p c (oo). 

A check of the functional form can be done by looking 
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FIG. 16: Pc(L) for the 24 NN model, using two definitions of 
the critical point: (1) the density at which P s = 1/2, and (2) 
the density at which the mean culling time peaks. 



at the width W of the transition, which we define to be 
the difference between the occupation probabilities that 
yield P s — 1/4 and P s = 3/4. The width as a function 
of L is shown in figure 17 for both wired and periodic 
boundary conditions. The same renormalization group 
formalism that implies | p c (L = oo) — p c (L) |oc L^ 1 ^ 
also implies W oc L~ x l v . However, Fig. 17 shows clear 
deviations from this power law form. This trend was 
starting to emerge in Ref. [11], though it was masked 
by the large error bars for the larger systems. In other 
words, the fit to the above assumption in Ref. [11] was 
premature. 
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FIG. 17: The widths for the 24NN model with both wired 
and periodic boundary condition, on a log-log plot. Both 
boundary conditions show clear deviations from power laws 
at large system sizes. 



Based on our above heuristic arguments for a con- 
nection between force-balance percolation and jamming 
percolation, we instead fit to the form for a diverging 
crossover length scale found by Toninelli, et al. [7, 8] 
(TBF) for the spiral model. Approaching the transi- 
tion from below, TBF identified a crossover length E, 
such that systems of size smaller than E are likely to 
have a spanning cluster, while those of size greater than 



S are exponentially unlikely to have a spanning clus- 
ter. TBF proved that E diverges at the transition faster 
than any power law. In Ref. [4], TBF argued that 
the upper and lower bounds scaled similarly, finding 
E - exp(-C(p -p c y), where fi = (1 - 0.64. 
(Both - = 0.63 and v\\ = 1.73 are from directed per- 
colation [46].) Based on the TBF formula, we propose 
that W oc (ln(L))- 1 /^. The widths are replotted with 
the appropriate scales for this fit in Figure 18, and while 
the curves still have deviations from perfect straight line 
behavior, the agreement is significantly better than for 
W oc L^ 1 ^ . From the width data, for wired bound- 
ary conditions, we extract fi — 0.39 ± 0.01; for periodic 
boundary conditions, fi = 0.45 ± 0.02. The stated error 
bar is purely statistical and does not take into account 
systematic effects that occur given the small range of sys- 
tem sizes, so one cannot necessarily rule out a link with 
overlapping directed percolation processes. We also note 
that in Ref. [4], TBF stated that their numerical data 
was consistent with fj, = 0.64 given the large systematic 
error bars. 
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FIG. 18: The widths for the 24NN model with both wired and 
periodic boundary condition, using a fitting form motivated 
by the TBF results. 



Now, if there is only one diverging lengthscale, then 
it can be argued that \p c (L) — p c (oo)\ should be propor- 
tional to W. So we choose the TBF fitting form and 
choose p c (oo) to minimize \ 2 f° r the best linear fit of 
ln(p c (oo) - Pc{L)) vs. ln(lnL), for L > 32. The op- 
timal fit, shown in figure 19, is obtained for p c {oo) = 
0.414 ± 0.008 for wired boundary conditions. (The error 
bar is obtained by testing which values of p c (oo) result 
in visibly noticeable curvature in the plot.) For periodic 
boundary conditions, we obtain p c (oo) = 0.425 ± 0.005. 
The straight line obtained in figure 19 does not make 
clear whether the TBF form is correct for our model, 
given that we have the freedom to choose p c (oo) in mak- 
ing the fit. From the fit with wired boundary conditions 
we obtain /i = 0.51 ± 0.09. Again, this result is to be 
compared with TBF's result of fj, = 0.64. 

From the two different results for fi, it appears that W 
and \p c (L = oo) — p c (L)\ do not scale in the same way. 
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There are two possible, but incompatible, explanations 
for the difference in the results for /i. One possibility is 
that the scaling of W and \p c (L = oo) — p c (L)\ is in fact 
the same in the thermodynamic limit, but that we obtain 
differing /i's for the finite system sizes studied due to the 
slow approach to the thermodynamic limit. This would 
be reminiscent of the extremely large finite-size effects 
seen in analogous systems, such as fc-core percolation, the 
Kob-Andersen model, and glassy dynamics. Numerically, 
this possibility is quite reasonable, given that the TBF 
formula for the growing length scale has an extra power 
law factor in front. A second possible explanation is that 
W and \p c (L — oo) —p c (L)\ in fact scale differently This 
would be possible if there are are two diverging length 
scales, rather than one; this would change the standard 
picture of a renormalization group fixed point controlled 
by a single parameter, so that W and \p c {L — oo) — p c (L)\ 
would no longer be forced to scale in the same manner. 
There are certainly two diverging lcngthscalcs in mean 
field fc-core percolation — one associated with those in the 
infinite cluster and one associated with sites that get re- 
moved in response to one random site being removed. 
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there is a jump in the order parameter at the transition, 
denoted by k c . In other words, in the thermodynamic 
limit, as soon as there is a spanning cluster, it occupies 
a finite fraction of the system. 

From Fig. 20, and the already-obtained result of p c — 
0.414 ± 0.008, we estimate the thermodynamic k c to be 
k c = 0.399 ± 0.011. Furthermore, for each individual 
curve it appears that just above the transition, k in- 
creases linearly with p, suggesting that the order param- 
eter exponent is j3 = 1. This result is to contrasted with 
the mean field fc-core percolation results and the repul- 
sive soft sphere simulations, where (3 = 1/2. 
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FIG. 20: Order parameter for the 24 NN model using wired 
boundary conditions. 
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FIG. 19: Fitting of p c (L) as a function of L using the TBF 
fitting form, for the 24 NN force-balance model with wired 
boundary conditions. We obtain the best fit with p c (oo) = 
0.414 ± 0.008. 



4- Order parameter 

We also investigate the order parameter, k, the frac- 
tion of sites in the infinite force-balance cluster. This 
is obtained for a given system size and initial occupa- 
tion density by generating configurations, keeping only 
those configurations in which the largest cluster is span- 
ning, and then finding the average fraction of sites in the 
largest cluster for those configurations. Plots are shown 
for L > 128 in figure 20. The curves lie on top of each 
other for large L, differing only in the minimum value of 
p needed for the numerical simulations to have a reason- 
able chance of obtaining occupied sites after culling. 

The minimum value of k is increasing, rather than de- 
creasing, with system size. We therefore assume that 



5. Correlation function 

We next look at correlation functions for this model, 
now using only periodic boundary conditions. For these 
simulations, for a given system of size L and density p, 
we generate states with exactly pi? occupied sites (rather 
than occupying sites independently). We only keep con- 
figurations that were nonempty after culling, and in each 
of these, keep only the largest cluster (this cluster was 
automatically spanning, since we use periodic boundary 
conditions). We then calculate the correlation function 
C{t) — the probability that sites a distance i apart are 
both occupied, minus an asymptotic constant, to be dis- 
cussed shortly. Moreover, we calculate correlations at all 
angles (i.e. not merely horizontal and vertical correla- 
tions) . 

With uncorrelated percolation, it is easy to quickly 
produce precise plots of the correlation function over sev- 
eral decades of distance, and confirm that the correlation 
function has a power law form at the critical point. How- 
ever, correlations in the force-balance model have signif- 
icantly greater sources of uncertainty. 

When we generate samples for uncorrelated percola- 
tion, each site occupation is independent, so the correla- 
tion function calculated for C(i) is largely (although not 
entirely) independent of that for C(J) with j ^ i. So 
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for uncorrelated percolation, when we look at different 
pairs of points in a single sample, we get somewhat inde- 
pendent estimates of the correlation function for a range 
of distances. That is, we quickly generate lots of inde- 
pendent data for the curve C(i). However, for correlated 
percolation models where there are no finite clusters, the 
results for C(i) and C(J), for j ^ i, are are highly corre- 
lated for a given sample, even for widely separated i and 
3- 

For a given amount of computer running time, the 
error bars in the correlation function are thus substan- 
tially larger than for uncorrelated percolation, and be- 
cause they are correlated, it is difficult to reliably extract 
the asymptotic value of C(i). This is important because 
to get C(i), we need to subtract off the probability that 
two distant sites (i — > oo) are occupied in the infinite 
system limit {L — > oo) after culling. For our models, we 
do not know the infinite system occupation probability 
after culling. (For uncorrelated percolation, it is trivially 
identical to the initial occupation probability, as there 
is no culling). The most straightforward solution is to 
subtract off the asymptotic value of C(i), since we know 
that the connected correlation function should approach 
as i — > oo. However, because the error bars at differ- 
ent i are correlated, there are also large error bars in the 
asymptotic value of C(i). The functional form of C(i) 
is highly sensitive to the constant subtracted off — over a 
limited distance range, a power law and an exponential 
can look very similar if the asymptotic value is shifted 
slightly. This means that that it is difficult to reliably 
determine the form of the correlation function. The cor- 
relation length, 



E^ 2 c(*) 



(i) 



is also very sensitive to a small shift in C(i). 

The end result is that it is difficult to get accurate re- 
sults for large system sizes, and we have limited our sim- 
ulations of correlation functions to L = 100, L = 200 and 
L — 400. The correlation length as a function of initial 
occupation density is shown for these three system sizes 
in figure 21. For densities p > 0.4, the correlations are 
clearly short-ranged, and are independent of the system 
size. For lower densities, the correlation length grows 
with L. The curve for L = 200 has a clear peak at 
p w 0.375; the L = 400 data appears to have a peak at 
a larger L, although the largish error bars make its loca- 
tion unclear. Thus, the data is somewhat suggestive of a 
lengthscale that grows with system size. 

The position of the peak gives yet another plausible 
measurement for p c (L). If we compare the p c (L) de- 
fined from the probability of spanning for L = 200, 
Pc(L) = 0.3671 ± 0.0007, which is lower than the esti- 
mated p c (L) from correlation length data. However, this 
discrepancy should vanish in the infinite system limit. 
Since the position of the peak is not so clear for L = 400 
it is more difficult to discern the trend. 



Because of the limited system sizes studied, it is 
somewhat difficult to determine with confidence whether 
the correlation functions have power law or exponential 
forms. The correlations are very short-ranged (£ < 5) 
for p > 0.4, indicating that the correlations are probably 
exponential (or a very steep power law) . At lower densi- 
ties, for most system sizes and densities, the correlation 
function is exponentially decaying at long distances; for 
example, see figure 22, showing the correlation functions 
for L = 400 and p = 0.37, 0.38, and 0.39. One excep- 
tion to the generally observed exponentialy behavior is 
te L = 200 correlation function for p = 0.395, which 
appears somewhat power-law-like, as shown in figure 23. 
The density of p = 0.395 differs slightly from the peak in 
the plot of £ for L = 200 in figure 21. For the L = 400 
system, the power-law-like trend of the correlation func- 
tion near the peak density is not as prominent and the 
correlation functions appear more exponential. 
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FIG. 21: Correlation lengths for the 24 NN force-balance 
model as a function of p, for L — 100 and L = 200. 
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FIG. 22: Correlation function for L = 400, for the 24 NN 
force-balance model, at p = 0.37, p = 0.38, and p — 0.39. For 
clarity, the p = 0.37 correlation function has been shifted up 
by 1.0, and the p — 0.39 correlation function has been shifted 
down by 1.0. 



These results should be compared with those of Parisi 
and Rizzo for k-core percolation in four dimensions [47]. 
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FIG. 23: Correlation function for L = 200, p = 0.395, for the 
24 NN force-balance model. The best fit line has a slope of 
-1.76 ±0.03. 



While they found that £ grew with decreasing p, they 
never obtained correlation lengths greater than 10. So 
they also found no power law correlations. They argue 
that the 4-core transition in four dimensions is an ordi- 
nary discontinuous transition with no diverging length- 
scales. We note that their systems had finite stable clus- 
ters and so may have different properties than the sys- 
tems considered in this paper. (The fact that their sys- 
tems had finite stable clusters also presumably means 
that they had smaller numerical correlations between 
C(i) and C(j) for | i — j | large, allowing them to deter- 
mine correlation functions for larger systems than here.) 



B. 16 NN Model and Spiral Model 

Now we present results for the 16 nearest neighbor 
force-balance model and spiral model. We begin with 
the mean culling time. The plot of the mean peak culling 
time as a function of L shows a similar 5/4 exponent as 
was measured for the 24NN model, a = 1.240± 0.025. In 
fact, all three two-dimensional models — the 24NN model, 
the 16NN model and the spiral model — are in quantita- 
tive agreement with another. We find a — 1.246 ± 0.027 
for the spiral model. These results suggest that similar 
processes underlie the culling in all three two-dimensional 
models. See Fig. 24. 

The probability of having a force-balance avalanche of 
size s, P(s), for the 16 NN model shows similar trends to 
the 24 NN model with a broad distribution of intermedi- 
ate sizes and a prominent peak for the largest sizes. See 
Fig. 25. The spiral model exhibits the same qualitative 
behaviour as well. 

Fig. 26 shows the probability of spanning for periodic 
boundary conditions. From this data we extract the 
width of the transition as a function of L, as depicted in 
Fig. 27. We plot both periodic and wired boundary con- 
ditions on a log-log scale to demonstrate the significant 
deivations from a power-law growing crossover length. 
We also plot the same data for the 24NN and spiral mod- 



els for comparison. As with the 24NN model, the fitting 
form is clearly not a power law in L. Figure 28 shows the 
widths using the fitting form motivated by the TBF re- 
sult. For the 16NN model, one can extract /i = 0.35±0.01 
for wired boundary conditions and /i = 0.42 ± 0.03 for 
periodic boundary conditions. 

In Fig. 29, similarly to Fig. 19, we choose p c (oo) so 
that p c (oo) — Pc(L) as a function of L is well fit by 
the TBF functional form. We obtain the best fit for 
p c (oo) = 0.497 ± 0.007 for wired boundary conditions, 
which gives /i = 0.47 ± 0.07. For periodic boundary con- 
ditions, p c (oo) = 0.502 ± 0.010 and fi = 0.70 ± 0.15. 
Once again, there is a discrepancy between the p ex- 
tracted from the width data and the /i extracted from 
the one-parameter fit. 

A similar fit (not shown) for the spiral model finds 
p c (oo) = 0.690 ± 0.008, and pL = 0.49 ± 0.08 for wired 
boundary conditions. For the spiral model, Toninclli, et 
al. [7, 8] have proven that p c (oo) is the same as for di- 
rect percolation, p c (oo) = 0.705, and found p = 0.64. 
Our numerical results for p c (oo) and fi are both within 
two error bars from their exact result, and a plot of 
ln(p c (oo) — p c (L)) versus ln(lnL) for the exact result of 
p c (oo) = 0.705 shows noticeable curvature. Also, the /x 
extracted from the width data for wired boundary condi- 
tions is n = 0.32 ± 0.01 and p = 0.38 ± 0.02 for periodic 
boundary conditions. 

Finally, Fig. 30 shows the order parameter K, clS cl func- 
tion of p for various system sizes for the 16NN model. 
As in the 24 NN model, the jump in n increases with 
increasing system size suggesting that the transition is 
discontinuous. It also appears that j3 = 1 just above the 
transition for each individual curve (though there is some 
overall curvature for the set of curves) . 
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FIG. 24: Peak of the mean culling time as a function of L for 
the all three models. For the 16NN model, the best fit has a 
slope of 1.240 ±0.025. For the spiral model, the best fit has a 
slope of 1.246 ± 0.027. Wired boundary conditions are used. 
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FIG. 25: Log-log plot of P(s) for the 16 NN model with L 
128 in the presence of periodic boundary conditions. 
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FIG. 27: The widths for all three models with both wired 
and periodic boundary condition, on a log- log plot. All curves 
show clear deviations from power laws at large system sizes. 
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FIG. 26: The probability of spanning for periodic boundary 
conditions for the 16 NN force-balance model. Error bars are 
too small to be seen. 
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FIG. 28: The widths for all three models with both wired and 
periodic boundary condition assuming the TBF form. 



C. Three-dimensional model 

Our tests of the 26NN three-dimensional model are 
not as extensive as in the two-dimensional case, particu- 
larly because we have no proof that p c < 1 in this case 
and our system size range is limited. From the prob- 
ability of spanning data using periodic boundary con- 
ditions, we have extracted W and p c (L). See Fig. 31. 
Deteriming p c (oo) with the TBF functional form yields 
p c (oo) = 0.433 ± 0.009 < 1 (not shown). Moreover, for 
the one-parameter p c (oo) fit, [i — 0.75 ± 0.14; from the 
width data, fi = 0.37 ± 0.01. In the three-dimensional 
case the discrepancy between the two values of // for 
the same boundary condition is greater than the two- 
dimensional force-balance models, where larger system 
sizes can be explored. For this three-dimensional model, 
one may expect a double logarthim as opposed a single 
logarithm with potentially /i = 0.52 using the values for 
three-dimensional directed percolation [46]. Given the 
small range of data, however, it is difficult to discrimi- 
nate between the two forms. 



Finally, the three-dimensional transition appears to be 
discontinuous in the context of the onset of the infinite 
force-balance cluster, as in the two-dimensional cases, 
since the jump in the order parameter increases with in- 
creasing system size. See Fig. 32. Moreover, simulations 
of the probability of having a force-balance avalanche size 
s shows the same trends as in Figs. 12 and 13. 

V. SUMMARY AND DISCUSSION 

While uncorrelated percolation is very well under- 
stood, models of correlated percolation are much less so. 
We have presented rigorous and numerical results on sev- 
eral force-balance percolation models to help narrow the 
gap. On the rigorous side, we have proven that p c < 1 for 
the two-dimensional models. This result places our inter- 
pretation of the numerical data for the two-dimensional 
models on sounder footing. A rigorous argument that 
the force-balance percolation transition is discontinuous, 
at least in two dimensions, is more difficult than in the 
case of jamming percolation. In the jamming percolation 
models the underlying mechanism driving the transition 
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FIG. 29: Fitting of p c (L) as a function of L using the TBF 
fitting form, for the 16 NN force-balance model with wired 
boundary conditions. We obtain the best fit with p c (oo) = 
0.497 ± 0.007. 



FIG. 31: Probability of having a spanning cluster in any of 
three directions, as a function of p, for different system sizes of 
length L, and periodic boundary conditions. This plot is for 
k — 4 with force-balance on the cubic lattice with 26 nearest 
neighbors. 
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FIG. 30: Order parameter for the 16 NN model using wired 
boundary conditions. 



is two disjoint, directed percolation processes. These two 
process scaffold upon one another such that the infinite 
cluster is bulky at the transition. The underlying mech- 
anism driving the transition in the force-balance models 
is presumably the same as argued in Sec. IIIB. Models of 
newly-constructed directed percolation processes will be 
key in constructing a rigorous argument for discontinuity 
in the force-balance case. 

Our numerical results for the force-balance avalanches 
and the onset of the infinite force-balance cluster for all 
three force-balance models strongly suggest a discontin- 
uous transition. For the force-balance avalanches, the 
probability of having an avalanche size s, is broad for in- 
termediate avalanche sizes. There is also a well-defined 
avalanche size at the tail of the distribution that becomes 
more prominent as the system size is increased and as p is 
decreased towards the transition. This suggests a bulky, 
discontinous, transition. Looking at the usual order pa- 
rameter also points towards a discontinuous transition in 
that the jump in the fractional size of the largest force- 
balance cluster at the transition increases with increasing 
system size. This trend was reported in Ref. [11] for the 
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FIG. 32: Order parameter plot for the three-dimensional 
model using periodic boundary conditions. 



24 NN model only. 

So now that we know there is a transition (at least 
for the two-dimensional models) and we surmise that the 
onset of the infinite force-balance cluster is discontinu- 
ous, just as in the case of jamming percolation, is there 
a quantitative connection with jamming percolation be- 
yond the heuristic arguments provided in Sec. IIIB? 

The results for the mean culling time quantitatively 
suggest that at least the two two-dimensional force- 
balance models and the spiral model are in the same uni- 
versality class. We obtained mean culling time exponents 
of a = 1.226±0.027, 1.240±0.025, and 1.246±0.027, for 
the 24 NN model, 16 NN model, and spiral model, respec- 
tively. Since all three exponents are within one standard 
deviation of 5/4, the equivalent exponent in the sandpilc 
model, there is potential for a sandpilc model-like RG 
treatment, at least for the culling dynamics for both sets 
of models. Further quantitative study of the distribution 
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of force-balance and spiral avalanche sizes will provide 
further insights. 

The crossover length data also indicates that for the 
system sizes studied, the mechanisms underlying jam- 
ming and force-balance percolation are the same for the 
two-dimensional models. The following table summa- 
rizes the values obtained for p c (oo) for the three two- 
dimensional models, and the associated values of \i: 



Model 


Boundary 




A* 


24 NN 


wired 


0.414 ±0.008 


0.51 ±0.09 


16 NN 


wired 


0.497 ±0.007 


0.47 ±0.07 


Spiral 


wired 


0.690 ± 0.008 


0.49 ±0.08 


24 NN 


periodic 


0.425 ± 0.005 


0.76 ±0.20 


16 NN 


periodic 


0.502 ±0.010 


0.70 ±0.15 


Spiral 


periodic 


0.701 ±0.011 


0.81 ±0.16 



We see that results for \i for all three models with wired 
boundary conditions are consistent with one another, al- 
though the error bars are large. So, while our measure- 
ments are not as precise, nor as accurate, as one would 
like (the former possibly indicated by the differing val- 
ues of [i obtained from the width data), the consistency 
between the three different models is readily apparent. 
In other words, the models are most likely in the same 
universality class. We have also included results for pe- 
riodic boundary conditions. For those, the results for 
p c (oo) are within the error bars of the results for wired 
boundary conditions; for periodic boundary conditions, 
the plots of ln(p c (oo) — p c (L)) versus ln(lnL) appear lin- 
ear for a wider range of p c (oo), resulting in significantly 
larger error bars. Finally, the p c 's are independent of the 
boundary condition, as expected. 

Our force-balance data is indeed inconsistent with a 
crossover length that diverges as a power law, and is 
just as consistent with the TBF fitting form as the spiral 
model data. With the present data, we are unable to 
conclude that the presumably two-dimensional model- 
independent value of n is independent of the bound- 
ary conditions, though with data for larger systems, 
and smaller error bars, such a trend should emerge. A 
stronger numerical test of the same underlying mecha- 
nism for jamming percolation and force-balance perco- 
lation would be to look for directed percolation using 
anisotropic finite-size scaling, as was done in Ref. [7]. 
We leave this for future work. 

Our correlation length data above the critical point for 
the 24 NN model suggests a correlation length that grows 
with system size, as opposed to a finite one. However, for 
the largest system sizes, we do not generally see the ex- 
pected power law correlations at the transition associated 
with this divergence. More work is needed to substanti- 
ate this trend, which would be inconsistent with a diverg- 
ing correlation length. Of course, a study of larger sys- 
tems might reveal a finite correlation length. Our current 
correlation length results are to be contrasted with 4-core 
percolation in four dimensions, where a finite correlation 



length of about 10 lattice spacings was found, suggest- 
ing a garden- variety type discontinuous transition driven 
by nucleation [47]. 4-core percolation in four dimensions 
contains finite clusters, which provide a backbone for nu- 
cleation. In force-balance percolation, however, there are 
no finite clusters so one may expect a more unusual dis- 
continuous transition. 

Therefore, the scenario for force-balance percolation 
that is most consistent with our data is that while the 
onset of the infinite force-balance cluster is discontinu- 
ous, there is an exponentially diverging crossover length- 
scale, and perhaps a diverging correlation lcngthscalc. 
Our limited data for the standard correlation length de- 
fined in the connected phase makes it difficult to discern 
any trend for growth, exponential or otherwise. In con- 
tinuous phase transitions, the crossover length and cor- 
relation length diverge in the same; way. With this more 
unusual transition, it is not necessarily obvious that the 
same behaviour should apply. Moreover, we find quanti- 
tative agreement with the dynamical exponent for sand- 
pile models not only for the force-balance models but for 
the spiral model as well, again, suggesting that they all 
are in the same universality class. Finally, we expect 
that three-dimensional versions of jamming percolation 
and force-balance percolation should exhibit similar be- 
haviour as well, as our data suggests. 

Though the usual order parameter, the fraction of sites 
in the infinite force-balance cluster, does not appear to 
be continuous, are there other candidates for an atypical 
continuous order parameter? A potential candidate is 
to look at the subset of the force-balance spanning clus- 
ter where the connectivity is marginal, i.e. 3-connected. 
However, we do not find any evidence for a fractal, span- 
ning 3-cluster at the force-balance percolation transition, 
nor for a fractal, spanning 4-cluster. Given that p c is 
around 0.4, each site has approximately 10 ncigbors at 
the transition. The 3-core condition is thus completely 
superseded by the vectorial constraint, at least for this 
lattice model with many nearest neighbors. While the 
dynamics of culling suggest a critical sandpilc-likc model 
for the removal of redundant sites, the removal of the 
marginal infinite cluster does not. So at this time, we 
have not discovered an order parameter which is contin- 
uous at the transition. 

While the jamming and force-balance percolation mod- 
els lead to a discontinuous transition (provably in the 
first case, and most likely in the latter), the fraction of 
the sites in the infinite cluster grows linearly in p — p c . 
It would certainly be interesting to uncover other finite- 
dimensional models that have a discontinuous transition 
in which the fraction of sites in the infinite cluster grows 
nonlinearly in p — p c above the transition; such mod- 
els would behave more like mean-field models. At this 
point, we know of no such models. It may be that finite- 
dimensional models of correlated connectivity percola- 
tion arc too simple to capture this aspect of jamming, 
and that one has to define forces on the network, as in 
rigidity percolation. We are currently working towards 
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this direction. 
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